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Abstract  -  A  method  is  presented  for  calculating  the  reconstruction  of  a  circularly  symmetric  two-dimensional 
function  from  its  projection,  a  relation  known  as  the  Abel  inversion.  This  technique  differs  from  techniques 
used  previously  by  using  integral  transforms  for  its  implementation.  The  frequency-domain  analysis  allows  for 
experimentally  obtained  data,  which  is  often  noisy  and  off-center,  to  be  dealt  with  in  a  systematic,  rational 
manner.  The  formulation  of  the  Abel  inversion  in  terms  of  transforms,  the  filtering  of  the  noise,  and  the  *«timite 
of  the  off-center  shift  are  discussed.  Sample  calculations  of  simulated  noisy  data  and  the  application  of  the  method 
to  an  image  of  a  laser  sustained  plasma  are  presented. 

L  Introduction 

Optical  emission  or  transmission  is  often  used  to  obtain  experimental  diagnostic  measurements  of  axially 
symmetric  sources  such  as  plasmas,  flames,  and  rocket  and  jet  engine  plumes.  The  line-of-sight  projection  of  the 
emission  or  transmission  coefficients  onto  a  one-dimensional  axis  can  be  used  to  determine  their  two-dimensional 
radial  structure.  As  is  the  case  in  most  inverse  problems,  however,  the  conditions  are  ill-posed,  and  so  considerable 
difficulty  can  arise  when  processing  real  experimentally  acquired  data  that  has  inherent  uncertainties. 

The  reconstruction  of  a  circularly  symmetric  two-dimensional  function  from  its  projection  onto  an  axis  is 
known  as  Abel  inversion  or  inverse  Abel  transformation  of  the  projection.  The  measured  intensity,  /(x),  is  given 
in  terms  of  emission  coefficients,  <(r),  through  the  Abel  transform1 

w 

* 

where  x  is  the  displacement  of  the  intensity  profile  and  r  is  the  radial  distance  in  the  source.  The  measured 
intensity  /(x)  is  the  one-dimensional  projection  of  the  two-dimensional  circularly  symmetric  function  having  c(r) 
as  a  radial  slice.  The  inversion  integral,  or  the  inverse  Abel  transform,  is  given  by, 


In  practice,  the  application  of  equation  (2)  is  made  difficult  because  of  the  singularity  in  the  integral  at  the 
lower  limit  and  because  the  derivative  of  the  projection  tends  to  greatly  enhance  the  noise  corrupting  the  data. 
Furthermore,  the  intensity  is  usually  not  available  as  a  continuous  function  and  so  is  known  only  at  discrete 
sample  points.  Several  approaches  based  upon  geometrical  techniques  or  numerical  methods  using  polynomial 
fits  have  been  used  to  perform  the  inversion.  Nestor  and  Olsen3  transformed  the  variables  according  to  r3  =  u 
and  x3  «  u  so  that  the  inversion  integral  can  be  approximated  by  a  simpler  sum.  Bockasten*  fitted  third  degree 
polynomials  to  the  data  points  and  approximated  the  integral  by  a  sum.  However,  these  methods  require  prior 
smoothing  of  the  data  and  are' not  considered  complete  in  themselves. 

Later,  least  squares  curve  fitting  methods  were  employed  and  were  found  to  yield  better  results  than  the  exact 
fit  methods  when  applied  to  noisy  data.  FVeeman  and  Kats4  used  a  single  polynomial  curve  to  fit  the  data.  Fourth 
order  polynomials  were  found  to  give  the  best  results  among  trials  using  up  to  twelfth  order  polynomials.  Cremers 
and  Birkebak*  compared  several  inversion  techniques  and  showed  that  the  least  squares  curve  fitting  techniques 
are  more  favorable  than  the  exact  fit  methods.  Short  descriptions  of  several  other  techniques  and  comparisons 
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of  these  various  techniques  can  be  found  in  reference  5.  Shelby6  divided  the  data  into  several  intervals  and 
used  a  least  squares  polynomial  fit  technique  in  each  interval  to  smooth  the  scattered  data.  The  inversion  was 
then  performed  analytically  and  summed  over  the  intervals  to  obtain  emission  coefficients.  Maldonado  et  a l.7 
expanded  e(r)  in  a  series  of  orthogonal  polynomials  and  derived  a  method  to  find  the  expansion  coefficients  from 
the  intensity  data. 

All  of  these  methods  have  drawbacks.  The  singularity  in  the  lower  limit  of  the  integral  causes  problems  for 
the  numerical  methods,  but  these  are  avoided  by  the  analytical  methods.  The  smoothing  techniques  used  are 
essentially  a  kind  of  lowpass  filtering  having  undetermined  filter  characteristics.  When  the  inversion  is  performed, 
the  spectral  characteristics  of  the  noise,  of  the  desired  signal  and  of  the  smoothing  algorithm  are  not  considered. 
Therefore,  the  possible  problems  of  loss  of  information  and  distortion  are  neglected.  Use  of  the  Abel  integral 
implies  an  assumption  that  the  input  data  will  be  symmetric,  but  the  measured  intensity  data  often  exhibit  some 
degree  of  asymmetry,  and  the  exact  axis  of  symmetry  is  not  usually  known  a  priori  Determination  of  the  axis  of 
symmetry  is  often  chosen  using  ad  hoe  methods.  The  smoothing  techniques  consume  a  large  amount  of  computer 
time  and  the  error  propagation  calculations  are  tedious6. 

We  have  developed  a  method  based  upon  integral  transforms  that  removes  many  of  the  difficulties  and 
uncertainties  associated  with  these  earlier  techniques  and,  by  using  the  fast  Fourier  transform  (FFT)  algorithm 
to  implement  the  procedure,  greatly  reduces  the  computation  time.  This  technique  is  substantially  different  from 
earlier  methods  in  its  use  of  transform  techniques  and  frequency  domain  analysis.  Several  principles  of  digital 
signal  processing  and  spectral  analysis  are  employed  to  solve  the  problems  associated  with  the  processing  of 
actual,  noise-corrupted,  asymmetric  data. 

In  Section  II  the  Abel  inversion  integral  is  reformulated  in  terms  of  the  Fourier  and  Hankel  transforms.  This 
reformulation  is  then  used  to  handle  many  of  the  ill-posed  conditions  encountered  when  analyzing  actual  noisy 
data,  and  the  methods  for  doing  this  are  outlined  in  Section  III.  Numerical  experiments  demonstrating  the  validity 
and  applicability  of  the  inversion  scheme  are  presented  in  Section  IV,  together  with  an  recent  application  of  the 
method  to  an  image  of  a  laser  sustained  plasma. 

II.  Reformulation  of  the  Abel  inversion 

The  reformulation  of  the  Abel  inversion  in  terms  of  integral  transforms  can  be  derived  beginning  with  equation 
(1).  With  the  substitution  r  =  y/x 2  +  y3,  the  projection  can  be  written 

/(*)  *  J  t{y/x* +y2)dy. 

-OO 

The  one-dimensional  Fourier  transform  of  equation  (3)  is 

oo  oo 

//  «( V*2  +  y7)exp(-j2xxq)dxdy. 

—  00  -00 

If  the  variables  of  integration  are  changed  from  cartesian  to  polar  coordinates,  it  can  be  shown  that* 


/T(/(z)]  =  2*j  rt(r)J0[2itrq)dr  (5) 

o 

where  J0(.)  denotes  the  sero-order  Bessel  function  of  the  first  kind.  The  right  hand  side  of  equation  (5)  is  the 
sero-order  Hankel  transform  of  «(r).  Thus,  the  emission  distribution  of  the  plasma  may  be  recovered  by  taking 
the  inverse  Hankel  transform  of  the  Fourier  transform  of  the  projected  intensity: 
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where  the  subscripts  R  and  I  correspond  to  the  real  and  imaginary  components  of  each  function.  Equation  (9b) 
allows  the  imaginary  portion  of  the  noise  spectrum  to  be  written  explicitly  in  terms  of  the  transform  of  the  data, 
D[k),  and  the  shift  of  the  data  from  the  origin,  Am. 

The  noise,  n(m),  was  further  assumed  to  be  an  uncorrelated  bandlimited  Gaussian  random  process.  Its 
spectral  components  could  thus  be  modeled  as  independent  Gaussian  random  variables*  each  with  probability 
density  function 

wWM)  ■  ^"'■1-^1'  (») 

Substituting  equation  (9b)  into  equation  (10)  and  taking  the  ensemble  probability  density  function  over  all  values 
of  Jb  as  the  product  of  the  individual  probability  density  functions  yields  the  conditional  probability  density  of 
jD(Jb)  given  Am  as 

1  w/3~* 

p(-P(*)|Am)  =  £  \Di{k)co»[2itk^mjM)  —  f?n(A:)stn(2jrfcAm/A/))2).  (11) 

Maximising  this  probability  density  function  with  respect  to  Am  provides  a  maximum  likelihood  estimator  for 
the  shift.  The  resulting  transcendental  equation  for  the  estimate,  Am,  is 

U/ 3-1 

53  k{Dn{k)Di(k)coa{4*kAm/M)  -  (/>£(*)  -  Z)?(fc)]«n(4*A:Am/M))  =  0.  (12) 
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Equation  (12)  was  solved  iteratively  to  determine  an  estimate  for  the  shift.  Each  spectral  component  of 
the  data  was  then  multiplied  by  the  phase  term  exp(— jlitk&m/M)  to  shift  the  data  back  to  its  proper  location 
centered  about  the  origin.  Following  this  multiplication,  the  imaginary  components  remaining  in  the  spectrum  of 
the  data  were  by  assumption  due  to  noise,  and  so  were  set  equal  to  tero  to  fully  symmetrize  the  data  and  further 
reduce  the  noise.  If  sufficient  data  is  available  to  provide  ensemble  averages,  this  noise  can  be  used  to  find  the 
statistical  properties  of  the  remaining  real  portion  of  the  noise  conrupting  the  signal,  but  this  approach  has  not 
been  taken. 

The  inversion  scheme  was  completed  by  taking  the  inverse  Hankel  transform  of  the  spectrum  of  the  filtered 
and  centered  data.  This  integral  transform  was  approximated  by  a  discrete  summation,  and  computation  time  was 
reduced  by  using  Candel's  algorithm10.  Further  reduction  in  computation  time  was  achieved  by  a  modification 
to  this  algorithm  which  requires  only  half  the  additions  of  the  Candel  method11. 

As  a  final  step  in  the  data  analysis,  a  study  of  the  effects  of  noise  on  the  errors  in  the  inverted  signal  was 
made.  If  the  input  noise  corrupting  the  intensity  data  is  assumed  to  be  a  wide-sense  stationary,  bandlimited 
random  process,  then  the  variance  of  the  inverted  output  noise  as  a  function  of  the  radial  coordinate  r  can  be 
expressed 

eo 

var[na[r)]  -  J (!r?)aG„(g)(J3(2srg)  +  Y*[2xrq)}dq  (13) 

o 

where  Gn[q)  is  the  power  spectral  density  of  the  input  noise.  For  the  often* assumed  case  of  bandlimited  white 
noise  with  variance  2  BN0,  and  bandwidth  B,  the  variance  of  the  output  noise  has  been  found  to  be  approximately 
equal  to 

vor(n0(r)j  »  (14) 

for  sufficiently  large  r.  These  equations  served  as  a  guide  in  the  design  of  filters  for  preprocessing  the  data,  and 
provided  a  benchmark  for  the  probable  errors  in  the  inverted  data. 


IV.  Numerical  and  experimental  results 


A  number  of  numerical  experiments  have  been  carried  out  using  both  noise-free  and  noisy  data  to  evaluate 
the  performance  of  the  integral  transform  method  implemented  using  discrete  transforms11.  These  numerical 
experiments  using  noise-free  Gaussian  and  semi-elliptic  data  established  the  accuracy  of  the  modified  Candel 
method  for  the  Hankel  transform.  Numerical  experiments  using  noisy  data  were  performed  to  establish  the 
validity  of  the  method  used  to  symmetrise  the  data  and  to  evaluate  the  effectiveness  of  the  optimal  filters  used 
to  smooth  the  data. 


A  noisy  trial  function  having  a  known  inversion  was  generated  by  sampling  a  Gaussian  function  defined  as, 


It  (m)  =  exp(-  ;  m  *  0, 1, ...,  127. 

The  Abel  inversion  of  I\  (m)  with  an  even  extension  in  negative  m  axis  is  theoretically  given  by,* 

€l^  *  I0^?e*p("l00)  ’  1  =  1 . 127' 


(15) 


(16) 


A  tero  mean  sequence  of  random  numbers  having  a  Gaussian  density  function  with  a  variance  of  .001  was 
generated  and  added  to  the  sequence  of  equation  (16)  which  was  then  shifted  by  two  integer  sample  points.  Also, 
the  sequence  was  sero  padded  so  that  the  filtering  by  discrete  Fourier  transforms  produces  a  linear  convolution 
filtering.  This  trial  sequence  is  shown  in  Figure  1. 


Figure  1.  Generated  noisy  and  shifted  trial  sequence  Ji(m). 


The  result  obtained  when  this  trial  sequence  was  inverted  without  filtering  or  symmetrizing  is  shown  in 
Figure  2.  The  same  trial  function  was  then  inverted  using  the  symmetrizing  procedure  and  an  optimal  filter 
having  a  passband  from  0.0  to  0.0525  and  a  stopband  from  0.0725  to  0.5.  The  result  from  this  inversion  is  shown 
in  Figure  3.  Other,  more  easily  generated  filters  were  found  to  produce  similarly  satisfactory  inversions. 

A  principal  advantage  of  this  integral  transform  method  is  its  computational  efficiency.  It  has  been  applied 
to  the  inversions  of  a  large  number  of  images  obtained  from  experiments  with  laser  sustained  plasmas13.  Each  of 
these  digitally  sampled  images  consisted  of  512  radial  scans  of  the  plasma  and  each  radial  scan  consisted  of  240 
sample  points.  The  inversion  of  512  sets  of  240-  point  data  consumed  about  eight  minutes  of  CPU  time  when 
processed  by  the  integral  transform  method,  in  contrast  to  a  20  hour  CPU  time  required  for  a  curve-fit  inversion 
using  the  method  described  in  reference  6.  The  data  analysis  was  carried  out  using  a  Masscomp  MC500  Micro 
Supercomputer.  The  symmetrizing  procedure  generally  consumed  about  five  to  seven  iterations  for  each  set  of 
data.  Contour  plots  of  the  intensity  and  the  emission  coefficient  for  one  of  these  images  are  given  in  Figure  4 
and  Figure  5.  This  particular  image  illustrates  the  performance  of  the  method  for  intensity  profiles  which  exhibit 
some  asymmetry  and  which  result  in  inversion  profiles  that  are  peaked  on-axis  (eg.  49  mm),  flat-topped  (eg.  45 
mm)  and  having  off-axis  peak s  (eg.  43  mm). 

IV.  Conclusions 

A  new  integral  transform  method  has  been  developed  for  the  Abel  inversion  of  measured  intensity  data  in 
the  presence  of  noise.  The  integral  transform  method  has  a  number  of  advantages  compared  with  the  existing 
curve-fit  techniques.  The  curve-fitting  methods  do  not  consider  the  spectral  characteristics  of  the  noise  and  the 
desired  information,  or  address  the  problem  of  symmetrizing  the  experimental  data.  The  Abel  inversion  or  the 
inverse  Abel  transform  is  performed  by  the  Fourier  transform  followed  by  the  inverse  Hankel  transform.  The 
efficiency  of  the  well  known  FFT  is  exploited  in  the  evaluation  of  both  the  Fourier  and  inverse  Hankel  transforms. 
This  fast  technique  is  made  applicable  to  noisy  and  off-center  experimental  data  by  frequency  domain  filtering 
and  a  maximum  likelihood  estimator  derived  from  the  assumed  noise  characteristics. 


*  o-o  o-o  o  computed 
-  exact 


technique. 
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Figure  4.  Contour  plot  of  the  projected  intensity  image  data  for  a  laser  sustained  plasma. 
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Figure  5.  Contour  plot  of  the  emission  coefficient  image  obtained  by  the  inversion  of  data  in  Figure  5. 
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